#########################################################################
# sample from the posterior distribution of Wakefield's baseline model
# for ecological inference with covariates R using linked C++ code in Scythe
# (King's Extended Model)
# 
# Incorporates parts of hierna.R and mcmcregress.r from MCMCpack ver x.xx
# Updated October 27, 2003
#
# Kevin Corder
#
#
##########################################################################

MCMChierEXC <- function(r0, r1, c0, c1, X, burnin=1000, mcmc=50000, thin=1,
                       nu0=1.0, delta0=0.5,
                       nu1=1.0, delta1=0.5,
                       beta0.start=NA, sigma20.start=NA,
                       beta1.start=NA, sigma21.start=NA,
                       b0=0, B0=0, b1=0, B1=0,
                       verbose=FALSE, tune=2.65316, seed=0, ...){

# Set up Regression quantities

   X <- as.matrix(X)         # X matrix
   N <- nrow(X)	             # number of observations      
   K <- ncol(X)              # number of covariates   
   P <- c0/(c0+c1)
   Y <- log(P/(1-P))
   

  # Error checking
  if (length(r0) != length(r1)){
    cat("ERROR: length(r0) != length(r1).\n")
    stop("Please check data and try again.\n")
  }

  if (length(r0) != length(c0)){
    cat("ERROR: length(r0) != length(c0).\n")
    stop("Please check data and try again.\n")
  }

  if (length(r0) != length(c1)){
    cat("ERROR: length(r0) != length(c1).\n")
    stop("Please check data and try again.\n")
  }

  if (length(r1) != length(c0)){
    cat("ERROR: length(r1) != length(c0).\n")
    stop("Please check data and try again.\n")
  }

  if (length(r1) != length(c1)){
    cat("ERROR: length(r1) != length(c1).\n")
    stop("Please check data and try again.\n")
  }

  if (length(c0) != length(c1)){
    cat("ERROR: length(c0) != length(c1).\n")
    stop("Please check data and try again.\n")
  }

  if (min((r0+r1) == (c0+c1))==0){
    cat("ERROR: rows and columns do not sum to same thing.\n")
    stop("Please check data and try again.\n")
  }

  if (burnin < 0 ){
    cat("ERROR: parameter burnin < 0.\n")
    stop("Please respecify and try again.\n")
  }

  if (mcmc < 0 ){
    cat("ERROR: parameter mcmc < 0.\n")
    stop("Please respecify and try again.\n")
  }

  if (mcmc%%thin != 0 ){
    cat("ERROR: parameter mcmc not evenly divisible by parameter thin.\n")
    stop("Please respecify and try again.\n")
  }

 if (nu0 <= 0 ){
    cat("ERROR: parameter nu0 <= 0.\n")
    stop("Please respecify and try again.\n")
  }

  if (nu1 <= 0 ){
    cat("ERROR: parameter nu1 <= 0.\n")
    stop("Please respecify and try again.\n")
  }

  if (delta0 <= 0 ){
    cat("ERROR: parameter delta0 <= 0.\n")
    stop("Please respecify and try again.\n")
  }

  if (delta1 <= 0 ){
    cat("ERROR: parameter delta1 <= 0.\n")
    stop("Please respecify and try again.\n")
  }

# starting values for beta0 error checking

   if (is.na(beta0.start)){ # use MLEs
      beta0.start <- matrix(coef(lm(P~X-1)), K, 1)
   }

   else if(is.null(dim(beta0.start))) {
      beta0.start <- beta0.start * matrix(1,K,1)  
      }
   else if((dim(beta0.start)[1] != K) || (dim(beta0.start)[2] != 1)) {
      cat("Starting value for beta0 not conformable.")
      stop("\n Please respecify and call the function again.\n")
      }

# prior for beta0 error checking

   if(is.null(dim(b0))) {
      b0 <- b0 * matrix(1,K,1)  
      }
   if((dim(b0)[1] != K) || (dim(b0)[2] != 1)) {
      cat("N(b0,B0^-1) prior b0 not conformable.")
      cat("\n Please respecify and call the function again.\n") 
      }  
 
   if(is.null(dim(B0))) {
      B0 <- B0 * diag(K)    
      }
   if((dim(B0)[1] != K) || (dim(B0)[2] != K)) {
      cat("N(b0,B0^-1) prior B0 not conformable.")
      stop("\n Please respecify and call the function again.\n")
      }  

# starting values for beta1 error checking

   if (is.na(beta1.start)){ # use MLEs
     beta1.start <- matrix(coef(lm(P~X-1)), K, 1)
   }

   else if(is.null(dim(beta1.start))) {
      beta1.start <- beta1.start * matrix(1,K,1)  
      }
   else if((dim(beta1.start)[1] != K) || (dim(beta1.start)[2] != 1)) {
      cat("Starting value for beta1 not conformable.")
      stop("\n Please respecify and call the function again.\n")
      }

# prior for beta1 error checking

   if(is.null(dim(b1))) {
      b1 <- b1 * matrix(1,K,1)  
      }
   if((dim(b1)[1] != K) || (dim(b1)[2] != 1)) {
      cat("N(b1,B1^-1) prior b1 not conformable.")
      cat("\n Please respecify and call the function() again.\n") 
      }  
 
   if(is.null(dim(B1))) {
      B1 <- B1 * diag(K)    
      }
   if((dim(B1)[1] != K) || (dim(B1)[2] != K)) {
      cat("N(b1,B1^-1) prior B1 not conformable.")
      stop("\n Please respecify and call the function() again.\n")
      }  
   
# sigma20 starting values error checking

   if (is.na(sigma20.start)){ # use MLE
     lm.out <- lm(P~X-1)
     sigma20.start <- var(residuals(lm.out))
   }

  else if(sigma20.start <= 0) {
      cat("Starting value for sigma20 negative.")
      stop("\n Please respecify and call the function again.\n")
      }   

# sigma21 starting values error checking

   if (is.na(sigma21.start)){ # use MLE
     lm.out <- lm(P~X-1)
     sigma21.start <- var(residuals(lm.out))
   }

  else if(sigma21.start <= 0) {
      cat("Starting value for sigma21 negative.")
      stop("\n Please respecify and call the function again.\n")
      }   
   
# prior for sigma20 error checking
 
  if(nu0 <= 0) {
      cat("IG(nu0/2,delta0/2) prior nu less than or equal to zero.")
      stop("\n Please respecify and call the function again.\n")
      }
   if(delta0 <= 0) {
      cat("IG(nu0/2,delta0/2) prior delta less than or equal to zero.")
      stop("\n Please respecify and call the function again.\n")      
      }  
 
# prior for sigma21 error checking
 
  if(nu1 <= 0) {
      cat("IG(nu1/2,delta1/2) prior nu less than or equal to zero.")
      stop("\n Please respecify and call the function again.\n")
      }
   if(delta1 <= 0) {
      cat("IG(nu1/2,delta1/2) prior delta less than or equal to zero.")
      stop("\n Please respecify and call the function again.\n")      
      }  

  # setup matrix to hold output from sampling
  ntables = length(r0)
  sample <- matrix(0, mcmc/thin, (ntables*2)+(K*2)+2)
  

  # call C++ code to do the sampling
  C.sample <- .C("hierEXC",
                 samdata = as.double(sample),
                 samrow = as.integer(nrow(sample)),
                 samcol = as.integer(ncol(sample)),
                 r0 = as.double(r0),
                 r1 = as.double(r1),
                 c0 = as.double(c0),
                 c1 = as.double(c1),
                 ntables = as.integer(ntables),
                 burnin = as.integer(burnin),
                 mcmc = as.integer(mcmc),
                 thin = as.integer(thin),
                 nu0 = as.double(nu0),
                 delta0 = as.double(delta0),
                 nu1 = as.double(nu1),
                 delta1 = as.double(delta1),
                 verbose = as.integer(verbose),
                 tune = as.double(tune),
                 seed = as.integer(seed),
                 accepts = as.integer(0),
                 Xdata = as.double(X),
	         Xrow = as.integer(nrow(X)),
	         Xcol = as.integer(ncol(X)),   
	         b0startdata = as.double(beta0.start),
                 b0startrow = as.integer(nrow(beta0.start)),
	         b0startcol = as.integer(ncol(beta0.start)),
	         b1startdata = as.double(beta1.start),
                 b1startrow = as.integer(nrow(beta1.start)),
	         b1startcol = as.integer(ncol(beta1.start)),
                 b0data = as.double(b0),
	         b0row = as.integer(nrow(b0)),
	         b0col = as.integer(ncol(b0)),   
	         B0data = as.double(B0),
	         B0row = as.integer(nrow(B0)),
	         B0col = as.integer(ncol(B0)),   
	         b1data = as.double(b1),
	         b1row = as.integer(nrow(b1)),
	         b1col = as.integer(ncol(b1)),   
	         B1data = as.double(B1),
	         B1row = as.integer(nrow(B1)),
	         B1col = as.integer(ncol(B1)),   
                 sigma20start = as.double(sigma20.start),
                 sigma21start = as.double(sigma21.start),
                 PACKAGE="MCMCpack"
                 )

  cat(" overall acceptance rate = ",
      C.sample$accepts / (C.sample$burnin+C.sample$mcmc) / ntables, "\n")

  sample <- matrix(C.sample$samdata, C.sample$samrow, C.sample$samcol,
                   byrow=TRUE)

  output <- mcmc2(data=sample, start=1, end=mcmc, thin=thin)
  p0names <- paste("p0table", 1:ntables, sep="")
  p1names <- paste("p1table", 1:ntables, sep="")
  b0names <- paste("beta0_",1:K,sep=" ")
  b1names <- paste("beta1_",1:K,sep=" ")
  varnames(output) <- c(p0names, p1names, b0names, "sigma^20", b1names, "sigma^21")

  attr(output, "title") <- "MCMCpack Wakefield's Hierarchical EI Model with covariates Posterior Density Sample"


  return(output)
  
}




















